Velocity Statistics in the Two-Dimensional Granular Turbulence 
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We studied the macroscopic statistical properties on the freely evolving quasi-elastic hard disk 
(granular) system by performing a large-scale (up to a few million particles) event-driven molecular 
dynamics systematically and found that remarkably analogous to an enstrophy cascade process in 
the decaying two-dimensional fluid turbulence. There are four typical stages in the freely evolving 
inelastic hard disk system, which are homogeneous, shearing (vortex), clustering and final state. 
In the shearing stage, the self-organized macroscopic coherent vortices become dominant. In the 
clustering stage, the energy spectra are close to the expectation of Kraichnan-Batchelor theory and 
the squared two-particle separation strictly obeys Richardson law. 
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The dynamics of granular materials becomes one of 
the most important topics in the studies of nonlinear, 
dissipative, and non-equilibrium statistical physics [1]. 
Granular media are collections of macroscopic particles 
with rough surfaces, dissipative, and frictional interac- 
tions. The granular systems require an energy source in 
order to be in a steady state and the external gravita- 
tional force much affect their dynamics. 

To focus on the dissipative features, smooth inelas- 
tic hard sphere (IHS) model are often used as an ideal 
model. The freely cooling granular fluid has been studied 
as an ideal dissipative particle system in the absence of 
external force. Since the system is only composed of an 
inelastic hard spheres and no relevant energy scale ex- 
ists, the restitution coefficient between collision particles 
is the only parameter in terms of non-equilibrium. The 
assumption of an inelastic hard sphere potential is also 
employed in kinetic theory, which facilitates comparisons 
between theory and simulation. In order to construct 
the theory of the macroscopic phenomenology in non- 
equilibrium dissipative particle system, an IHS model is 
the most promising as a microscopic model. A linear sta- 
bility analysis of hydrodynamic equations for IHS model 
has revealed that the initial spatially homogeneous cool- 
ing state is unstable to the formation of vortices and clus- 
ters. The shearing (vortex) and cluster instabilities were 
theoretically predicted and were tested by molecular dy- 
namics (MD) simulations [2]. 

Since the total energy is monotonically decreasing in 
the freely evolving process, a steady state in terms of en- 
ergy fluctuation can be realized by scaling the velocity of 
the entire particle to the total energy remaining constant 
[3]. Here, we introduce the new-scaled time tg which is 
described by 



dt 

WV 



Pit) = vm/m, 



(0.1) 



where T(t) is the average kinetic energy per particle as 
a function of usual time t. The great advantage of this 
scaling operation is that the trajectories of particles don't 
change compared with the non-scaled case in case of hard 
sphere system. Therefore, one simply replaces the usual 
time t with the new-scaled time ts. This operation is 
the same as the well-known velocity scaling method [4] 
in the usual MD simulation, in which the velocities of all 
particle Vi(t) are scaled following each collision by the 
factor P{t) (i.e. P{t)-Vi{t)) and the total energy is kept 
fixed strictly all over the time. 

The 2D turbulence in nature is a large-scale fluid mo- 
tion in the atmosphere or ocean dynamics on earth. The 
most remarkable feature of 2D turbulence is described by 
the enstrophy cascade dynamics, which is completely dif- 
ferent from that of 3D turbulence represented by the K41 
theory. The existence of enstrophy cascade process was 
originally proposed by Kraichnan [5] and Batchelor [6]. 
The theory expected that the enstrophy injected at a pre- 
scribed scale is dissipated at smaller scales, undergoing a 
cascading process at a constant enstrophy transfer rate; 
this led to predicting a spectrum for the energy, in a 
range of scales extending from the injection to the dissi- 
pative scale. The granular kinetic energy spectrum in the 
connection with the fluid turbulence was firstly pointed 
out by Taguchi [7] in 2D granular vibrated beds. He ob- 
tained the results of the k^^^^ spectrum in his simulation 
with a few hundred particles. Another important nature 
of 2D turbulence is the self-organized coherent vortices, 
which develop into larger ones through the merging pro- 
cess between vortices with same sign of circulation. 

In this letter, we especially focus on the velocity statis- 
tics and statistical laws on the fluid turbulence. To spec- 
ify what is the universal characters in dissipative system 
both macroscopic equation and microscopic dissipative 
particle, we performed extensive event-driven molecular 
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dynamics simulation systematically on a freely cooling 
process in 2D IHS model with velocity scaling thermo- 
stat. We found a strong similarity between 2D IHS model 
and 2D NS fluid turbulence. 

The freely 2D IHS (granular) model is so simple that 
the system is completely characterized with only three 
dimensionless parameters: the restitution coefficient r, 
the total number of disks A'', and the packing fraction 
p. The system size in the unit of disk diameter d is 
L/d = ^J-kN jvj^. All the disks are identical, namely 
the system is monodispersed. The collision is instanta- 
neous and only binary collisions occur. When two disks, 
i and j, with respective velocities Vj and Vj collide, the 
velocities after the collision, v- and v^-, are given by 

v^=Vi-^(l + r)[n.(v,-v,)]n (0.2) 
v^- = V,- + ^(1 + r)[n • (v, - v,)]n, (0.3) 

where n is the unit vector parallel to the relative posi- 
tion of the two colliding disks in contact. Our system 
consists of more than 250 thousands hard disks (up to a 
few million) placed in a square box with periodic bound- 
aries without no external force. To perform such a large- 
scale simulation, we implemented the simple and efficient 
event-driven algorithm, which can actually simulate more 
than a few million particles even in the personal com- 
puter [8]. Initially, the system is prepared as the equi- 
librium state by the long enough preliminary run with 
the restitution coefficient r = 1, in which the density is 
uniform and the disk velocities are Maxwell-Boltzmann 
distribution. The packing fraction and the restitution 
coefficient (i^, r) were varied from dilute to dense and 
within shearing regime, which is estimated by the crite- 
rion of McNamara and Young [9] , respectively. In case of 
(iV, v) = (1024, 0.25), McNamara and Young have found 
the final states have three typical states, which are ki- 
netic, shearing, and collapse regime. However, the spatio- 
temporal structure of shearing regime have not known yet 
especially at the macroscopic level. 

The criterion of kinetic-shearing boundary in Ref. [9] 
is based on the results of Jenkins and Richman [10], 
in which the high wave number cutoff for the unstable 
shear modes were derived. On the contrary, the shearing- 
collapse boundary is estimated by ID theoretical analy- 
ses for inelastic collapse, which are known as the phe- 
nomenon on the number of collisions during a finite time 
diverges. By using the the regime criterion described 
by McNamara and Young [9], one can find both regime 
boundaries become close to the unity in the thermody- 
namic limit. These theoretical expectation indicate that 
the system always unstable even in the quasi-elastic limit. 
This fact implies the important conjecture discussed later 
when we consider a large-scale IHS model as the macro- 
scopic fluid model. 

In the large-scale simulations, the restitution coeffi- 
cients within shearing regime become quasi-elastic (r ~ 



1). All our simulations by changing various parameters 
within the shearing regime, the system evolves to the 
final steady state after several stages. As described in 
McNamara and Young [9], we can calculate the packing 
fraction J^(x), velocity u(x) = {ux{x),Uy{x)), and tem- 
perature T(x) at any point x. Figure 1 represents typ- 
ical evolving process for four normalized properties as a 
function of new-scaled time tg in 2D IHS model. During 
relaxation process, four stages can be distinguished. Af- 
ter the homogeneous cooling state (HCS) continues for a 
certain time from initial thermal equilibrium state(first 
stage), the short-range velocity correlation for each time 
(i.e. both pre- and postcoUisional velocity correlation) 
within the distance u(x) • u(x') sharply deviate from 
zero (solid line in Fig.l), where x' is location around x 
with a distance of disk diameter d (second stage) . Note 
that there are several works on the existence for the short- 
range velocity correlations even in HCS [11,12]. Our sim- 
ulation also seems to show that velocity correlation 'grad- 
ually' increase from the beginning of simulation in first 
stage. Therefore, it might to be difficult to determine the 
exact time between first and second stage. In this stage, 
coherent vortices sclf-organizc and the coherent vortices 
develop into larger ones through the mutual confiuence 
and the merging process among vortices with same sign of 
circulation (Fig. 2 (a)). These self-organized vortices are 
found firstly by McWilliams [13] in the direct numerical 
simulation (DNS) of 2D NS fiuid turbulence. In the third 
stage, the density fiuctuation Vrms (dotted line in Fig.l), 
which is calculated the square root of the space average 
(//(x) — ly)'^, exhibit instability compressive flow (Fig. 2 
(b)). Finally, steady state is realized (fourth stage). In 
the final steady states, the spatial correlation of inelastic 
hard disks, which gradually increases from the beginning 
of simulation, might be reached beyond the system size 
and begin to interfere with each other though the peri- 
odic boundary condition. Our simulations, in the shear- 
ing regime, there is no sign for the inelastic hard disks as- 
sembling to one cluster during the simulation time. This 
is because the shear mode expanded to the whole system 
might be stable though the periodic boundary condition. 
Therefore, we call them 'final steady states'. We found 
there are four characteristic final steady states, which 
are shear (laminar, oscillating, and turbulent) and vor- 
tex (one pair of vortices with opposite sign of circiflation) 
fiows, by changing both packing fraction and the resti- 
tution coefficient within shearing regime, systematically. 
Vortex flows of final pattern are also observed in the DNS 
simulation for 2D incompressible turbulence with the pe- 
riodic boundary condition. 

The velocity anisotropy A = X](mx(x)^ — 
Uy(x)^)/(?ij;(x)^ +Uy{x)'^), lu wffich One can distinguish 
the final state as a shear (^ = —1 or 1) or a vortex 
{A = 0) fiows, and the enstrophy Z = J2 l'^(x)P) where 
a;(x) = rot u(x), are also plotted by dot dashed line 
and dashed line in Fig. 1, respectively. We confirmed 
that the total vorticity J2x^i^) throughout the 

simulation. 
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The 2D fluid turbulence has different characters on the 
statistical law between forced and freely decaying case. 
However, in the granular case, no systematic considera- 
tion seems to exist yet. In the previous studies related 
to energy spectrum in granular material, energy injec- 
tion (thermostat) arc driven by the vibration cycle [7] 
and the periodical-stochastic thermostat [14]. These en- 
ergy injections resemble in the studies of forced fluid tur- 
bulence. On the other hand, velocity scaling thermo- 
stat [3], in which the system is driven continuously, is 
thought as corresponding to a freely decaying case, be- 
cause the statistics itself don't change by introducing the 
new-scaled time tg. Actually, we found the energy spec- 
tra (power spectra of velocity field) in 2D IHS model with 
the velocity scaling thermostat are close to the expecta- 
tion of Kraichnan-Batchelor theory {E{k) ~ k~^) after 
the third stage (Fig. 3). Therefore, our simulations show 
the enstrophy cascade, as is expected by the theory for 
freely decaying 2D fluid turbulence. We have confirmed 
this power law by changing several different system size. 
In Fig. 3, we can estimate the characteristic spatial scale 
{kd ~ 0.3) for minimal dissipative domain (such as Kol- 
mogorov scale in the fluid turbulence), which is composed 
of about a thousand disks. 

The first quantitative phenomenological observation in 
developed turbulence was shown by Richardson [15], in 
which the two fluid particle separation R= jr^ — [ obeys 
power law ((i?^) ~ t^). We also show that the time (ts) 
dependence of the space-averaged squared two-disk sep- 
aration in 2D IHS model strictly obeys Richardson dis- 
persion law in the third stage (Fig. 4). In the inset of 
Fig. 4, the enstrophy evolution is also plotted in terms 
of the new-scaled time tg. The enstrophy seems to de- 
cay power-law behavior in the second stage, in which the 
coherent vortices self-organize. However, since the sec- 
ond stage itself is relatively short, this behavior needs 
further conflrmation. As the intermittency of vorticity, 
Mc Williams found that the probability distribution func- 
tion of vorticity significantly deviate from Gaussian [13]. 
By calculating flatness of vorticity {fa, = (w^)/(w^)^) in 
2D IHS model, our simulations also show the deviation 
from Gaussian after third stage. 

How should we understand the obtained results? The 
different points between fluid turbulence and granular 
turbulence are compressibility, the origin of dissipation 
(i.e. viscosity and inelasticity between particles) and the 
ratio of particle size and system size. Is the granular tur- 
bulence close to fluid turbulence in the thermodynamic 
limit? Let us assume the extreme condition that is dense, 
thermodynamic and clastic limit. In this situation, a lit- 
tle amount of dissipation in the system always results in 
unstable state. We also obtained the fact that the veloc- 
ity correlation length (Kolmogorov scale) becomes larger 
in dense (i.e., quasi- incompressible) system, but less than 
system size when we consider the thermodynamic limit. 
Therefore, this extreme condition seems to really corre- 
spond to NS fluid turbulence. 

In this letter, we showed the remarkable similar as- 



pects on the statistical law of vorticity between 2D IHS 
(granular) turbulence and 2D NS fluid turbulence. These 
results were obtained by only solving a simple Newton's 
equation system for inelastic hard disks in terms of event- 
driven scheme. From the microscopic dynamics of inelas- 
tic hard disk to the macroscopic fluid, it is important to 
study the origins of statistical law for turbulence at the 
microscopic level, but there are very few studies so far 
from this point of view. The discussion for three lim- 
its (dense, thermodynamic, and elastic) might make a 
connection between 2D granular turbulence and 2D fluid 
turbulence. 
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FIG. 1. The time evolution of various properties from a 
equilibrium, which are density fluctuation (dotted line), veloc- 
ity correlation within a short distance (solid lino), anisotropy 
of total velocity (dot dashed line), and enstrophy (dashed 
line), respectively. The 0(ts) indicated four normalized prop- 
erties as a function of t,. The parameters are fixed at 
(r,iV,!/) = (0. 99109, 512^,0.60) during the simulation. An in- 
set in the below- left-hand corner shows the early stage of time 
evolution when shearing and clustering instability appears. 



FIG. 3. Energy spectra of the veloc- 
ity field are plotted for each stage. The parameters are at 
[r, N, i^)=(0.99109, 512^, 0.60). 

FIG. 4. The evolving squared two-particle separation in 
terms of new-scaled time ts is plotted. An inset in the up- 
per-left hand corner shows the time dependence of enstrophy 
decay The parameters are at (r, A^, i/)=(0.99109, 512^ 0.60). 



FIG. 2. The typical snapshots for coherent vortices and 
turbulent clustering patterns are shown. (a) The abso- 
lute vorticity field after 1200 collisions per particle with 
(r, AT, i^)=(0.99452, 640000, 0.65). The self-organized coher- 
ent vortices in the vorticity field grow spontaneously from 
the initial equilibrium state, (b) The density field after 
3590 collisions per particle in inelastic hard disk system with 
(r,Ar,j/)=(0.99725, 2560000, 0.65). The turbulent compres- 
sive flow appears in density field. 
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